《Gibbs Sampling for the UniniTiated》阅读笔记结构：
这篇是中篇，介绍一个非常简单的朴素贝叶斯文档模型生成的例子，用来说明Gibbs Sampler具体是如何构造的。
文档生成的建模过程
首先我们有一批文档，文档里面有很多单词，这些单词都是无顺序可交换的（词袋模型），这些文档分成两类，类标签为0或者1。给予一篇未标记的文档\(W_j\)，我们要做的工作就是预测文档的类标签是\(L_j=0\)还是\(L_j=1\)。为了方便起见，我们定了类标签所表示的类\(\mathbb{C}_0={W_j|L_j=0}\)和\(\mathbb{C}_1={W_j|L_j=1}\)。一般来说预测这种事都是选择最有可能发生的，即找到\(W_j\)的后验概率\(P(L_j|W_j)\)最大的标签\(L_j\)。使用贝叶斯公式 \[\begin{equation} \begin{split} L_j=\arg \max \limits_{L}P(L|W_j)& =\arg \max \limits_{L}\frac{P(W_j|L)P(L)}{P(W_j)}\\& =\arg \max \limits_{L} P(W_j|L)P(L) \\\end{split} \end{equation}\] 因为分母\(P(W_j)\)与\(L\)无关所以删去了。 通过贝叶斯公式的转换，我们可以想象这些文档的生成过程。首先，我们选择文档的类标签\(L_j\);假设这个过程是通过投硬币完成的（正面概率为\(\pi=P(L_j=1)\) )，正式地来说，就是服从贝努利分布 \[\begin{equation}L_j \sim Bernoulli(\pi)\end{equation}\] 然后，对于文档上\(R_j\)个“词位”中的每一个，我们根据一个概率分布\(\theta\)，随机独立地抽样一个词\(w_i\)。因为每个类生成词的\(\theta\)分布都不同，所以应该有\(\theta_1\)和\(\theta_2\)，具体地生成词的时候，我们根据文档的标签\(L_j\)来决定由哪个类来生成 \[\begin{equation} W_j \sim Multinomial(R_j,\theta_{L_j}) \end{equation}\]
先验
上面提到的参数\(\pi\)和\(\theta\)还有各自的先验知识。我们假设参数\(\pi\)是从一个参数为\(\gamma_{\pi_1}\)和\(\gamma_{\pi_0}\)的Beta分布中采样出来的。这里\(\gamma_\pi=<\gamma_{\pi_1},\gamma_{\pi_0}>\) 被称作超参数，因为它们是先验模型的参数，是为了确定模型的参数而存在的。类似地，就像Beta分布是贝努利分布的共轭先验一样，\(\theta\) 参数的确定来自于参数为\(\gamma_\theta\)的Dirichlet分布。公式如下： \[\begin{equation}\pi \sim Beta(\gamma_\pi)\end{equation}\]\[\begin{equation}\theta \sim Dirichlet(\gamma_\theta)\end{equation}\] 选择Beta分布和Dirichlet分布作为先验分布，是为了数学上的方便。所以的参数定义见图1和图2(这个概率图非常好看，是从别人的博客对这篇文章的介绍中借来的，来自于[3]）.
状态空间以及初始化
状态空间。在前面提到过一样，Gibbs 抽样器要做的是遍历模型中\(k\)维状态空间，对里面所有的变量采样。在当前模型的状态空间中存在以下变量：
1个标量变量\(\pi\)
2个向量变量\(\theta_0\)和\(\theta_1\)。
N篇文档所对应的二元标签向量变量\(\mathbf{L}\)。
还有每篇文档里的词向量\(W_j\)，但是他们是被观察到的数据，值已经知道了（这就是为啥图2中的概率图\(W_{jk}\)是实心的）。
初始化。初始化就是一篇文档的生成过程，不细讲了。
产生联合分布
对于之前提到的每次循环\(t=1\ldots T\)，我们都像上篇的式子(11)中提到的那样，利用其他变量的条件分布对状态空间中的每个变量更新。具体过程如下：
写出所有变量的联合分布
简化联合分布的公式
确定上篇的式（11）中条件分布的公式
给出采样器伪码的最终形式
写出联合分布
根据模型，我们可以得到整个文档的联合分布\(P(\mathbb{C},\mathbf{L},\pi,\theta_0,\theta_1;\gamma_{\pi_1},\gamma_{\pi_0},\gamma_{\theta})\)。 根据模型的生成过程和概率图，我们可以将联合分布分解: \[\begin{equation} P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})P(\mathbf{L}|\pi)P(\theta_0|\gamma_{\theta})P(\theta_1|\gamma_{\theta})P(\mathbb{C}_0|\theta_0,\mathbf{L})P(\mathbb{C}_1|\theta_1,\mathbf{L}) \end{equation}\] 然后将一个一个地来解释这些式子(从引用了一部分）：
\(P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})\)。这个式子是根据超参数为\(\gamma_{\pi_1}\)和\(\gamma_{\pi_0}\)的Beta分布中采样出\(\pi\)的公式，根据Beta 分布的定义，概率为 \[\begin{equation}P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})=\frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})}\pi^{\gamma_{\pi_1}-1}(1-\pi)^{\gamma_{\pi_0}-1} \end{equation}\] 因为右边式子的\(\frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})}\) 是个常量(设为\(c\))，而整个联合分布最终要把所有的常量归一化，所以在这里先不管它，把式子写作 \[\begin{equation} P(\pi|\gamma_{\pi_1},\gamma_{\pi_0}) \propto \pi^{\gamma_{\pi_1}-1}(1-\pi)^{\gamma_{\pi_0}-1} \end{equation}\]
\(P(\mathbf{L}|\pi)\)。第二个式子是由参数\(\pi\)或\(1-\pi\)生成标签变量1或0的过程，要注意一共有\(N\)篇文档，每篇文档的标签生成都是相互独立的： \[\begin{equation} \begin{split}P(\mathbf{L}|\pi)& =\prod^N_{n=1}\pi^{L_n}(1-\pi)^{N-L_n} \\ &=\pi^{C_1}(1-\pi)^{C_0} \\ \end{split} \end{equation}\] 其中\(C_0\)是和\(C_1\)分别是标签为0和1的文档数目。
\(P(\theta_0|\gamma_{\theta})\)和\(P(\theta_1|\gamma_{\theta})\)。这个式子是根据超参数为\(\gamma_{\theta}\)的Dirichlet分布中采样出\(\theta_0\)和\(\theta_1\)的词分布的过程。因为\(\theta_0\)和\(\theta_1\)相互独立，生成的过程是一样的，所以为了简单起见，我们暂时舍去类的下标，公式如下： \[\begin{equation} \begin{split} P(\theta|\gamma_{\theta})&=\frac{\Gamma{(\sum^V_{i=1}\gamma_{\theta_i})}}{\prod^V_{i=1}\Gamma{(\gamma_{\theta_i})}}\prod^V_{i=1}\theta_i^{\gamma_{\theta_i}-1}\\& = c'\prod^V_{i=1}\theta_i^{\gamma_{\theta_i}-1} \\& \propto \prod^V_{i=1}\theta_i^{\gamma_{\theta_i}-1} \\\end{split} \end{equation}\] \(\gamma_{\theta_i}\)代表了向量\(\gamma_{\theta}\)第\(i\)维的值。类似地，\(\theta_i\)代表向量\(\theta\)地\(i\)维的值，其物理意义就是该分布上第\(i\)个词被生成的概率。\(c'\)是另一个归一化因子。
\(P(\mathbb{C}_0|\theta_0,\mathbf{L})\)和\(P(\mathbb{C}_1|\theta_1,\mathbf{L})\)。这俩式子是具体地生成该类的所有文档中的所有词的过程。首先要求来看对于单独一个文档\(n\)，产生所有word也就是\(W_n\)的概率。假设对于某个文档，\(\theta=(0.2,0.5,0.3)\)，意思就是word1产生概率为0.2，word2产生概率为0.5，word3的概率为0.3。假如这个文档里word1有2个，word2有3个，word3有2个，则这个文档的产生概率就是(0.20.2)(0.50.50.5)(0.30.3)。所以按照这个道理，一个文档生成的联合概率如下： \[\begin{equation}P(W_n|\mathbf{L},\theta_{L_n})=\prod^V_{i=1}\theta_i^{W_{ni}} \end{equation}\] 这里\(\theta_{L_n}\)代表了文档\(n\)所在的类（\(\theta_0\)或\(\theta_1\)），\(\theta_i\)像之前提过的那样，代表产生第\(i\)个词的概率，而指数\(W_{ni}\)为该词出现的频率。现在我们有一篇文档的生成概率了，然后我们把一个类别下面的所有文档都合并起来(因为都是相互独立的么) \[\begin{equation}\begin{split} P(\mathbb{C}_x|\mathbf{L},\theta_x)& =\prod_{n\in \mathbb{C}_x} \prod^V_{i=1}\theta_i^{W_{ni}} \\ &=\prod^V_{i=1}\theta_{x,i}^{\mathcal{N}_{\mathbb{C}_x(i)}}\end{split}\end{equation}\] 其中\(x\)代表了类的标号（这里是0或1），\(\mathcal{N}_{\mathbb{C}_x}\)代表了类标号为\(x\)的所有文档生成词\(i\)的数目。
先验选择和简化联合分布
然后我们将先验知识 \(P(\mathbf{L}|\pi)\)与观察到的evidence\(P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})\) 相乘，我们可以看到估计的参数\(\pi\) 是如何被观察到的数据所影响的。利用式子(7) 和(9) 我们得到 \[\begin{equation} \begin{split} P(\pi|\mathbf{L};\gamma_{\pi_1},\gamma_{\pi_0}) &=P(\mathbf{L}|\pi)P(\pi|\gamma_{\pi_1},\gamma_{\pi_0}) \\& \propto [\pi^{C_1}(1-\pi)^{C_0}][\pi^{\gamma_{\pi_1}-1}(1-\pi)^{\gamma_{\pi_0}-1}] \\& \propto \pi^{C_1+\gamma_{\pi_1}-1}(1-\pi)^{C_0+\gamma_{\pi_0}-1} \\ \end{split} \end{equation}\] 对于\(\theta\),类似地我们结合(10）和（12）,有 \[\begin{equation} \begin{split}P(\theta|W_n;\gamma_\theta)& =P(W_n|\theta)P(\theta|\gamma_\theta) \\& \propto \prod^V_{i=1}\theta_i^{W_{ni}} \prod^V_{i=1}\theta_i^{\gamma_{\theta_i}-1} \\& \propto \prod^V_{i=1}\theta_i^{W_{ni}+\gamma_{\theta_i}-1} \\\end{split}\end{equation}\] 然后从一篇文章推广到某类下的所有文章，我们有 \[\begin{equation} P(\theta_x|\mathbb{C}_x;\gamma_\theta) \propto \prod^V_{i=1}\theta_{x,i}^{\mathcal{N}_{\mathbb{C}_x}(i)+\gamma_{\theta_i}-1} \end{equation}\] 最后，我们把所有的因式都乘起来，再令\(\mu=<\gamma_{\pi_1},\gamma_{\pi_0},\gamma_{\theta}>\),写出简化后的全联合概率: \[\begin{equation} P(\mathbb{C},\mathbf{L},\pi,\theta_0,\theta_1;\mu) \propto \pi^{C_1+\gamma_{\pi_1}-1}(1-\pi)^{C_0+\gamma_{\pi_0}-1} \prod^V_{i=1}\theta_{0,i}^{\mathcal{N}_{\mathbb{C}_0}+\gamma_{\theta_i}-1} \theta_{1,i}^{\mathcal{N}_{\mathbb{C}_1}+\gamma_{\theta_i}-1} \end{equation}\]
将隐含变量\(\pi\)积出
为了进一步地简化模型，我们可以对参数\(\pi\)求积分，从而消去这个参数。 \[\begin{equation} \begin{split}& P(\mathbb{C},\mathbf{L},\theta_0,\theta_1;\mu)\\ =&\int_\pi P(\mathbb{C},\mathbf{L},\pi,\theta_0,\theta_1;\mu) d\pi \\=&\int_\pi P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})P(\mathbf{L}|\pi)P(\theta_0|\gamma_{\theta})P(\theta_1|\gamma_{\theta})P(\mathbb{C}_0|\theta_0,\mathbf{L})P(\mathbb{C}_1|\theta_1,\mathbf{L}) d\pi \\ =&P(\theta_0|\gamma_{\theta})P(\theta_1|\gamma_{\theta})P(\mathbb{C}_0|\theta_0,\mathbf{L})P(\mathbb{C}_1|\theta_1,\mathbf{L}) \int_\pi P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})P(\mathbf{L}|\pi) d\pi \\\end{split}\end{equation}\] 这里我们对积分式子内的\(\pi\)求积分 \[\begin{equation} \begin{split}&\int_\pi P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})P(\mathbf{L}|\pi) d\pi \\ =&\frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})}\pi^{\gamma_{\pi_1}-1}(1-\pi)^{\gamma_{\pi_0}-1} \pi^{C_1}(1-\pi)^{C_0} d\pi \\=& \frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})} \int_\pi \pi^{C_1+\gamma_{\pi_1}-1}(1-\pi)^{C_0+\gamma_{\pi_0}-1} d\pi \\\end{split}\end{equation}\] 要注意这里式子右边的积分部分其实是尚未归一化的参数为\(C_1+\gamma_{\pi_1}\)和\(C_0+\gamma_{\pi_0}\) 的Beta分布求积分，所以我们只需要先补上Beta(\(C_1+\gamma_{\pi_1}\),\(C_0+\gamma_{\pi_0}\))归一化的常数项即可: \[\begin{equation}\frac{\Gamma(C_0+C_1+\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(C_1+\gamma_{\pi_1})\Gamma(C_0+\gamma_{\pi_0})}\end{equation}\] 令\(N=C_0+C_1\)，我们得到（注意这里应为Beta归一化的常数的倒数！！！原文中写反了） \[\begin{equation} \int_\pi P(\pi|\gamma_{\pi_1},\gamma_{\pi_0})P(\mathbf{L}|\pi) d\pi=\frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})} \frac{\Gamma(C_1+\gamma_{\pi_1})\Gamma(C_0+\gamma_{\pi_0})}{\Gamma(C_0+C_1+\gamma_{\pi_1}+\gamma_{\pi_0})} \end{equation}\] 最后得出 \[\begin{equation}P(\mathbb{C},\mathbf{L},\theta_0,\theta_1;\mu) \propto \frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})} \frac{\Gamma(C_1+\gamma_{\pi_1})\Gamma(C_0+\gamma_{\pi_0})}{\Gamma(C_0+C_1+\gamma_{\pi_1}+\gamma_{\pi_0})}\prod^V_{i=1}\theta_{0,i}^{\mathcal{N}_{\mathbb{C}_0}(i)+\gamma_{\theta_i}-1} \theta_{1,i}^{\mathcal{N}_{\mathbb{C}_1}(i)+\gamma_{\theta_i}-1}\end{equation}\] 这里很容易想到的是，既然我们可以通过对\(\pi\)求积分来简化联合分布，为什么不用这种方法把里面的参数全积分了呢？详细请见下篇的解释。
构建Gibbs Sampler
根据Gibbs 抽样的定义，我们每次都从条件分布中抽样出\(Z_i\)的新值： \[\begin{equation} P(Z_i|z_1^{(t+1)},\ldots,z^{(t+1)}_{i-1},z^{(t)}_{i+1},\ldots,z^{t}_k) \end{equation}\] 替换成具体的例子，在分配\(L_1^{(t+1)}\)的值时，我们需要计算条件分布 \[\begin{equation}P(L_1|L_2^{(t)},\ldots,L_N^{(t)},\mathbb{C},\theta_0^{(t)},\theta_1^{(t)};\mu)\end{equation}\] 然后我们算\(L_2{(t+1)}\)的值， \[\begin{equation}P(L_2|L_1^{(t+1)},L_3^{(t)},\ldots,L_N^{(t)},\mathbb{C},\theta_0^{(t)},\theta_1^{(t)};\mu)\end{equation}\] 以此类推直到计算完\(L_N^{(t+1)}\)为止。然后我们计算\(\theta_0\)的值, \[\begin{equation}P(\theta_0^{(t+1)}|L_1^{(t+1)},L_2^{(t+1)},\ldots,L_N^{(t+1)},\mathbb{C},\theta_1^{(t)};\mu)\end{equation}\] 最后是\(\theta_1\) \[\begin{equation}P(\theta_1^{(t+1)}|L_1^{(t+1)},L_2^{(t+1)},\ldots,L_N^{(t+1)},\mathbb{C},\theta_0^{(t)};\mu)\end{equation}\] 在循环\(t\)开始的时候，我们拥有一些信息，这些信息包括：每个文档中词的数目，标签为0的文档数，标签为1的文档数，所有标签为0的文档的词数目，所有标签为1的文档的词数目，每个文档的当前标签，当前的\(\theta_1\)和\(\theta_0\)值，等等。当我们想得到文档\(j\)的新标签时，我们暂时地移除所有当前文档的信息（包括词数目和标签信息），然后通过余下的信息得出\(L_j=0\) 的条件概率和\(L_j=1\)的条件概率，最后根据这俩概率的相对比例采样得到新的\(L_j{(t+1)}\)。对\(\theta\)采样时也是如此。
对文档标签\(\mathbf{L}\)采样
接下来我们具体地看看那如何对文档标签采样。根据之前的条件概率公式，我们得到 \[\begin{equation} \begin{split} P(L_j|\mathbf{L}^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)& =\frac{P(L_j,W_j,\mathbf{L}^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)}{P(\mathbf{L}^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)} \\ &=\frac{P(\mathbf{L},\mathbb{C},\theta_0,\theta_1;\mu)}{P(\mathbf{L}^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)} \end{split} \end{equation}\] \(\mathbf{L}^{(-j)}\)是除了\(L_j\)外所有的文档标签，然后\(\mathbb{C}^{(-j)}\)是除了\(W_j\)外所有的文档集合。这个分布有2个结果，\(L_j=0\)或\(L_j=1\)。 要注意这里分子就是式(32)中的全联合概率分布。而分母仅仅除去了文档\(W_j\) 的信息。然后我们来看看除去了该文档对整个式子造成了什么影响。
式（32）中第一个因式\(\frac{\Gamma(\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(\gamma_{\pi_1})\Gamma(\gamma_{\pi_0})}\)是个常数，与\(W_j\)无关，分子分母都有一个，所以计算的时候被约掉了。第二个式子是 \[\begin{equation} \frac{\Gamma(N+\gamma_{\pi_1}+\gamma_{\pi_0})}{\Gamma(C_1+\gamma_{\pi_1})\Gamma(C_0+\gamma_{\pi_0})} \end{equation}\] 现在我们来看看去掉了一个文档\(W_j\)，发生了什么变化。首先，语料的大小从\(N\)减小到\(N-1\)。 如果文档\(L_j=0\)，那么\(C_0^{(-j)}=C_0-1,C_1^{(-j)}=C_1\)。反之如果文档\(L_j=1\)，那么\(C_1^{(-j)}=C_1-1,C_0^{(-j)}=C_0\)。我们令\(x \in {0,1}\)这样就可以统一以上的情况，\(C_x^{(-j)}=C_x-1\)。然后我们重写分子和分母的这两个式子，从 \[\begin{equation} \frac{\frac{\Gamma(C_1+\gamma_{\pi_1}}{\Gamma(N+\gamma_{\pi_1}+\gamma_{\pi_0}))\Gamma(C_0+\gamma_{\pi_0})}}{\frac{\Gamma(C_1^{(-j)}+\gamma_{\pi_1})\Gamma(C_0^{(-j)}+\gamma_{\pi_0})}{\Gamma(N+\gamma_{\pi_1}+\gamma_{\pi_0}-1)}} \end{equation}\] 到 \[\begin{equation} \frac{\Gamma(C_x+\gamma_{\pi_x})\Gamma(N+\gamma_{\pi_1}+\gamma_{\pi_0}-1)}{\Gamma(N+\gamma_{\pi_1}+\gamma_{\pi_0})\Gamma(C_x+\gamma_{\pi_x}-1)} \end{equation}\] 利用伽马函数\(\Gamma(a+1)=a\Gamma(a)\)的性质，我们简化上面的式子最终得到 \[\begin{equation} \frac{C_x+\gamma_{\pi_x}-1}{N+\gamma_{\pi_1}+\gamma_{\pi_0}-1} \end{equation}\] （注意！！！这里原文又写错了，分子项少减了一个1）这样我们就把伽马函数消去了。再来看剩下的式子 \[\begin{equation} \prod^V_{i=1}\theta_{0,i}^{\mathcal{N}_{\mathbb{C}_0}(i)+\gamma_{\theta_i}-1} \theta_{1,i}^{\mathcal{N}_{\mathbb{C}_1}(i)+\gamma_{\theta_i}-1} \end{equation}\] 当我们去掉了文档\(W_j\)之后，可以发现文档不属于的另一类是不会有变化的，会上下约掉。而文档属于的那一类，\(C_x^{(-j)}=C_x-1\)，与文档\(W_j\)无关的都被約掉，只剩下\(W_j\)中的词，最后 \[\begin{equation} \prod^V_{i=1} \frac{\theta_{x,i}^{\mathcal{N}_{\mathbb{C}_x}(i)+\gamma_{\theta_i}-1}}{\theta_{x,i}^{\mathcal{N}_{\mathbb{C}_x^{(-j)}}(i)+\gamma_{\theta_i}-1}}=\prod^V_{i=1} \theta_{x,i}^{W_{ji}} \end{equation}\] 然后把（42）式和（44）式相乘，就得到了最后的公式，对于\(x\in{0,1}\) \[\begin{equation} P(L_j=x|\mathbf{L})^{(-j)},\mathbb{C}^{(-j)},\theta_0,\theta_1;\mu)=\frac{C_x+\gamma_{\pi_x}-1}{N+\gamma_{\pi_1}+\gamma_{\pi_0}-1}\prod^V_{i=1} \theta_{x,i}^{W_{ji}} \end{equation}\] 这个式子清楚地表现了标签是如何被选择的。这个式子中，前半部分其实只有\(C_x\)是变量，所以如果\(C_0\)大，则算出来的\(P(L_j=0)\)的概率就会大一点，所以下一次Lj的值就会倾向于C0，反之就会倾向于C1。 而后半部分，是在判断当前\(\theta\)参数的情况下，这些词\(W_j\)是否符合该参数所对应的分布（表现为似然值更大），然后确定整个文档是更倾向于C0还是C1。
最后（终于结束了！），我们从得到的条件概率中抽样：
令\(value0\)为式(34)\(x=0\)时的值。
令\(value1\)为式(34)\(x=1\)时的值。
令概率分布为\(<\frac{value0}{value0+value1},\frac{value1}{value0+value1}>\)。
根据这个分布进行贝努利实验（投两边概率不同的硬币）得到新的\(L_j^{(t+1)}\)的值。
对\(\theta\)采样
接下來看看如何对\(\theta_0\)和\(\theta_1\)采样。因为这俩参数都是相互独立的，我们再一次删去他们的下标，统一地看做\(\theta\)。将式(25)中与当前\(\theta\)无关的参数删去，我们可以观察到 \[\begin{equation}P(\theta|\mathbb{C},\mathbf{L});\mu) \propto P(\mathbb{C},\mathbf{L}|\theta)P(\theta|\mu)\end{equation}\] 其实这时候式子里只剩下\(\theta_{x,i}^{\mathcal{N}_{\mathbb{C}_x}(i)+\gamma_{\theta_i}-1}\)而已。可以看到，因为先验分布\(P(\theta|\mu)\)是共轭的Dirichlet分布，所以后验分布也是Dirichlet分布，只是每个词项\(i\)的参数变成了\(\mathcal{N}_{\mathbb{C}_x}(i)+\gamma_{\theta_i}\)而已，我们定义\(V\) 维的向量\(\mathbf{t}\),令\(t_i=\mathcal{N}_{\mathbb{C}_x}(i)+\gamma_{\theta_i}\)，就有 \[\begin{equation}\theta \sim Dirichlet(\mathbf{t})\end{equation}\] 那么如何从一个Dirichlet分布采样出\(\theta\)呢？（见的11.1.2的Reject Sampling）为了从参数为\(<\alpha_1,\ldots,\alpha_V>\)的\(V\)维Dirichlet分布中抽样出一个随机向量\(\mathbf{a}=<a_1,\ldots,a_V>\)，我们只要从\(V\)个伽马分布中各自采样一个独立的样本\(y_1,\ldots,y_V\)，其中每个伽马分布的密度函数为 \[\begin{equation}Gamma(\alpha_i,1)=\frac{y_i^{\alpha_i-1}e^{-y_i}}{\Gamma(\alpha_i)}\end{equation}\] 然后令\(a_i=y_i/\sum^V_{j=1}y_j\)即可。
利用已有标签的文档
我们可以利用已经被标签过的文档来，只需要不对这些文档采样新标签即可。这些被标注的文档会作为背景信息为算法服务。
整个过程
最终的算法图：
需要注意的是每次新的标签\(L_j\)产生的时候，都会对接下来的文档标签数目统计产生影响，这就是Gibbs Sampler的本质。
从Gibbs sampler产生值
Gibbs Sampling在每个循环都会产生变量的值，像之前提到过的那样，在理论上，变量\(Z_i\)可以通过\(T\)次获得的值近似 \[\begin{equation}\frac{1}{T}\sum^T_{t=1}z^{(t)}_i\end{equation}\] 但是实际中一般不直接这样用。
收敛和burn-in迭代
根据选择初始值的不同，Gibbs sampler需要一定的迭代次数才能保证点\(<z_1^{t},z_2^{t},\ldots,z_k^{t}>\)都是从马尔科夫链的平稳分布中生成的（换句话说就是马尔科夫链需要一定的次数才能收敛）。为了避免在这之前的估计对结果产生的影响，一般都丢弃\(t<B\)之前的结果，之前的阶段就被称为“burn-in”阶段，所以取平均值的时候是从\(B+1\)次到\(T\)次的。
自相关和lag
式（38）里的近似假设\(Z_i\)的那些样本都是相互独立的，而事实上我们知道他们不是，因为新的点都是从前面的点所给的条件所生成的。这个问题被称作自相关(autocorrelation)。为了避免这个问题，很多Gibbs Sampling在实现的时候取每\(L\)个值的平均值，这个\(L\)被称作滞后（lag，我也不知道怎么翻译）。具体的探讨请见原文。
还有多链问题和超参数的选择问题，以及原文推荐的一些有用的文章，这里就不详述了。
参考文献：
- 《Pattern Recognition and Machine Learning》
- 《LDA数学八卦》
- http://www.xperseverance.net/blogs/2013/03/1682/
- 《Gibbs Sampling for the UniniTiated》